
## @knitr MODELS_SETUP_OVERALL_MID

rm(list=ls())

# from \url{https://www.r-bloggers.com/identifying-the-os-from-r/}
get_os <- function(){
  sysinf <- Sys.info()
  if (!is.null(sysinf)){
    os <- sysinf['sysname']
    if (os == 'Darwin')
      os <- "osx"
  } else { ## mystery machine
    os <- .Platform$OS.type
    if (grepl("^darwin", R.version$os))
      os <- "osx"
    if (grepl("linux-gnu", R.version$os))
      os <- "linux"
  }
  tolower(os)
}

get_machine  <- function(){
    sysinf  <- Sys.info()
    if (!is.null(sysinf)){
        my.machine  <- sysinf['nodename']
    }
    tolower(my.machine)
}


if(get_os()=="windows" & get_machine()=="mompracen") {
    path <- setwd('C:/Users/giaco/Documents/MEGAsync/')
} else if (get_os()=="windows" & get_machine()!="mompracen") {
    path <- setwd('C:/Users/grieco/Dropbox/')
} else if (get_machine()=="aramis") {
    path <- setwd('/home/giacomo68/Dropbox/')
} else {
    path <- setwd('/home/giacomo/Dropbox/')
}


library(survival)
library(xtable)
library(ggplot2)
library(english)
library(Publish)
library(showtext)
library(extrafont)
library(BFpack)
library(mice)
library(mitools)
library(sandwich)
library(lmtest)


Data <- read.csv('gc-jg-project/kampala/data/data-kampala.csv',
                 stringsAsFactors=FALSE)


Data$country_id <- as.numeric(factor(Data$country_name))

m <- 10

imp <- mice(Data[,c("kampala",
                 "sidea.1991",
                 "sideb.1991",
                 "my.mil.cap.lag",
                 "v2x_polyarchy.lag",
                 "v2x_libdem.lag",     
                 "v2x_partipdem.lag",
                 "v2x_delibdem.lag", 
                 "v2x_egaldem.lag",
                 "gdp.cap.lag",
                 "pop.lag",
                 "common.law",
                 "year",
                 "rule.law.lag",
                 "aid.recipient",
                 "aid.recipient.q",
                 "All_NGO.lag",
                 "country_id")], m = m, maxit = 10, seed = 1234)

models <- list()
model_results <- list()
merged_imputations <- list()
m1.model_results <- list()
m2.model_results <- list()
m3.model_results <- list()
m4.model_results <- list()
m5.model_results <- list()
m6.model_results <- list()

BF.m1.mid <- list()
BF.m2.mid <- list()
BF.m3.mid <- list()
BF.m4.mid <- list()
BF.m5.mid <- list()
BF.m6.mid <- list()

BF.mid.table <- list()
BF.mid.bayes.factors <- list()

for (i in 1:m) {
    complete_data <- complete(imp, i)

    complete_data$.imp <- i
    
    merged_data <- merge(complete_data, 
                         Data[,c("country_id","year","country_name",
                                 "kampala.t0","kampala.t1",
                                 "sidea.1945",
                                 "sideb.1945",
                                 "sidea.force.1945",
                                 "sideb.force.1945",
                                 "sidea.no.force.1945",
                                 "sideb.no.force.1945",
                                 "sidea.force.1991",
                                 "sideb.force.1991",
                                 "sidea.no.force.1991",
                                 "sideb.no.force.1991",
                                 "mid.1991",
                                 "mid.1945",
                                 "mid.force.1945",
                                 "mid.no.force.1945",
                                 "mid.force.1991",
                                 "mid.no.force.1991"
                                 )], 
                        by = c("country_id", "year"),  # Note: quoted variable names
                        all.x = TRUE)  # Keep all rows from complete_data
    
    merged_imputations[[i]] <- merged_data


    current_data <- merged_imputations[[i]]

    m1 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
              mid.1991
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log(gdp.cap.lag)
            + log(pop.lag)
            + common.law
            + cluster(country_id)
           ,data=current_data)

     m1.vcov_cluster <- vcovCL(m1, cluster = current_data$country_id)

     m1.model_results[[i]] <- list(
        coefficients = coef(m1),
        variance = m1.vcov_cluster,
        n = nrow(current_data)
       )
    

    m2 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                mid.1945
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log(gdp.cap.lag)
            + log(pop.lag)
            + common.law
            + cluster(country_name)
           ,data=current_data)

    m2.vcov_cluster <- vcovCL(m2, cluster = current_data$country_id)

    m2.model_results[[i]] <- list(
        coefficients = coef(m2),
        variance = m2.vcov_cluster,
        n = nrow(current_data)
       )
    


    m3 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    mid.force.1991
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log(gdp.cap.lag)
                + log(pop.lag)
                + common.law
                + cluster(country_name)
               ,data=current_data)

    m3.vcov_cluster <- vcovCL(m3, cluster = current_data$country_id)

    m3.model_results[[i]] <- list(
        coefficients = coef(m3),
        variance = m3.vcov_cluster,
        n = nrow(current_data)
       )
    
    

    m4 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    mid.force.1945
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log(gdp.cap.lag)
            + log(pop.lag)
            + common.law
            + cluster(country_name)
           ,data=current_data)

    m4.vcov_cluster <- vcovCL(m4, cluster = current_data$country_id)

    m4.model_results[[i]] <- list(
        coefficients = coef(m4),
        variance = m4.vcov_cluster,
        n = nrow(current_data)
       )
    


    m5 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    mid.no.force.1991
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log(gdp.cap.lag)
                + log(pop.lag)
                + common.law
                + cluster(country_name)
               ,data=current_data)
    
    m5.vcov_cluster <- vcovCL(m5, cluster = current_data$country_id)

     m5.model_results[[i]] <- list(
        coefficients = coef(m5),
        variance = m5.vcov_cluster,
        n = nrow(current_data)
       )
    


    m6 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    mid.no.force.1945
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log(gdp.cap.lag)
                + log(pop.lag)
                + common.law
                + cluster(country_name)
               ,data=current_data)
    
    m6.vcov_cluster <- vcovCL(m6, cluster = current_data$country_id)

    m6.model_results[[i]] <- list(
        coefficients = coef(m6),
        variance = m6.vcov_cluster,
        n = nrow(current_data)
       )
    


    BF.m1.mid[[i]] <- BF(m1,hypothesis="mid.1991<0;mid.1991=0;mid.1991>0",complement=TRUE)
    BF.m2.mid[[i]] <- BF(m2,hypothesis="mid.1945<0;mid.1945=0;mid.1945>0",complement=TRUE)
    BF.m3.mid[[i]] <- BF(m3,hypothesis="mid.force.1991<0;mid.force.1991=0;mid.force.1991>0",complement=TRUE)
    BF.m4.mid[[i]] <- BF(m4,hypothesis="mid.force.1945<0;mid.force.1945=0;mid.force.1945>0",complement=TRUE)
    BF.m5.mid[[i]] <- BF(m5,hypothesis="mid.no.force.1991<0;mid.no.force.1991=0;mid.no.force.1991>0",complement=TRUE)
    BF.m6.mid[[i]]  <- BF(m6,hypothesis="mid.no.force.1945<0;mid.no.force.1945=0;mid.no.force.1945>0",complement=TRUE)
    
    BF.mid.table[[i]] <- rbind(BF.m1.mid[[i]]$PHP_confirmatory,
                               BF.m2.mid[[i]]$PHP_confirmatory,
                               BF.m3.mid[[i]]$PHP_confirmatory,
                               BF.m4.mid[[i]]$PHP_confirmatory,
                               BF.m5.mid[[i]]$PHP_confirmatory,
                               BF.m6.mid[[i]]$PHP_confirmatory)
    
    BF.mid.bayes.factors[[i]] <- rbind(BF.m1.mid[[i]]$BFmatrix_confirmatory,
                                       BF.m2.mid[[i]]$BFmatrix_confirmatory,
                                       BF.m3.mid[[i]]$BFmatrix_confirmatory,
                                       BF.m4.mid[[i]]$BFmatrix_confirmatory,
                                       BF.m5.mid[[i]]$BFmatrix_confirmatory,
                                       BF.m6.mid[[i]]$BFmatrix_confirmatory)

    
}


m1.pooled_results <- MIcombine(
    results = lapply(m1.model_results, function(x) x$coefficients),
    variances = lapply(m1.model_results, function(x) x$variance)
)

m2.pooled_results <- MIcombine(
    results = lapply(m2.model_results, function(x) x$coefficients),
    variances = lapply(m2.model_results, function(x) x$variance)
)

m3.pooled_results <- MIcombine(
    results = lapply(m3.model_results, function(x) x$coefficients),
    variances = lapply(m3.model_results, function(x) x$variance)
)

m4.pooled_results <- MIcombine(
    results = lapply(m4.model_results, function(x) x$coefficients),
    variances = lapply(m4.model_results, function(x) x$variance)
)

m5.pooled_results <- MIcombine(
    results = lapply(m5.model_results, function(x) x$coefficients),
    variances = lapply(m5.model_results, function(x) x$variance)
)

m6.pooled_results <- MIcombine(
    results = lapply(m6.model_results, function(x) x$coefficients),
    variances = lapply(m6.model_results, function(x) x$variance)
)

# View final results
summary(m1.pooled_results)
summary(m2.pooled_results)
summary(m3.pooled_results)
summary(m4.pooled_results)
summary(m5.pooled_results)
summary(m6.pooled_results)


coef.m1 <- m1.pooled_results$coefficients
se.m1 <- sqrt(diag(m1.pooled_results$variance))
z.m1 <- coef.m1/se.m1
p.m1 <- 2*pnorm(-abs(z.m1))
conf.m1 <- confint(m1.pooled_results)
m.m1 <- cbind(coef.m1,se.m1,p.m1,conf.m1)

coef.m2 <- m2.pooled_results$coefficients
se.m2 <- sqrt(diag(m2.pooled_results$variance))
z.m2 <- coef.m2/se.m2
p.m2 <- 2*pnorm(-abs(z.m2))
conf.m2 <- confint(m2.pooled_results)
m.m2 <- cbind(coef.m2,se.m2,p.m2,conf.m2)

coef.m3 <- m3.pooled_results$coefficients
se.m3 <- sqrt(diag(m3.pooled_results$variance))
z.m3 <- coef.m3/se.m3
p.m3 <- 2*pnorm(-abs(z.m3))
conf.m3 <- confint(m3.pooled_results)
m.m3 <- cbind(coef.m3,se.m3,p.m3,conf.m3)

coef.m4 <- m4.pooled_results$coefficients
se.m4 <- sqrt(diag(m4.pooled_results$variance))
z.m4 <- coef.m4/se.m4
p.m4 <- 2*pnorm(-abs(z.m4))
conf.m4 <- confint(m4.pooled_results)
m.m4 <- cbind(coef.m4,se.m4,p.m4,conf.m4)

coef.m5 <- m5.pooled_results$coefficients
se.m5 <- sqrt(diag(m5.pooled_results$variance))
z.m5 <- coef.m5/se.m5
p.m5 <- 2*pnorm(-abs(z.m5))
conf.m5 <- confint(m5.pooled_results)
m.m5 <- cbind(coef.m5,se.m5,p.m5,conf.m5)

coef.m6 <- m6.pooled_results$coefficients
se.m6 <- sqrt(diag(m6.pooled_results$variance))
z.m6 <- coef.m6/se.m6
p.m6 <- 2*pnorm(-abs(z.m6))
conf.m6 <- confint(m6.pooled_results)
m.m6 <- cbind(coef.m6,se.m6,p.m6,conf.m6)

                 

rownames(m.m1) <- rownames(m.m2) <- rownames(m.m3) <-
    rownames(m.m4) <- rownames(m.m5) <- rownames(m.m6) <- 
    c("MID",
      "Milex per cap",
      "Polyarchy",
      "GDP cap (log)",
      "Population (log)",
      "Common law")

colnames(m.m1) <- colnames(m.m2) <- colnames(m.m3) <-
    colnames(m.m4) <- colnames(m.m5) <- colnames(m.m6) <-
    c("b","std.err","p","X2.5.","X97.5.")

t.m1 <- data.frame(round(m.m1,3))
t.m1$Variables <- rownames(m.m1)
rownames(t.m1) <- NULL
t.m1$type <- "Overall,\nfrom 1991"

t.m2 <- data.frame(round(m.m2,3))
t.m2$Variables <- rownames(m.m2)
rownames(t.m2) <- NULL
t.m2$type <- "Overall,\nfrom 1945"

t.m3 <- data.frame(round(m.m3,3))
t.m3$Variables <- rownames(m.m3)
rownames(t.m3) <- NULL
t.m3$type <- "Force,\nfrom 1991"

t.m4 <- data.frame(round(m.m4,3))
t.m4$Variables <- rownames(m.m4)
rownames(t.m4) <- NULL
t.m4$type <- "Force,\nfrom 1945"

t.m5 <- data.frame(round(m.m5,3))
t.m5$Variables <- rownames(m.m5)
rownames(t.m5) <- NULL
t.m5$type <- "No force,\nfrom 1991"

t.m6 <- data.frame(round(m.m6,3))
t.m6$Variables <- rownames(m.m6)
rownames(t.m6) <- NULL
t.m6$type <- "No force,\nfrom 1945"

t.combined <- rbind(t.m1,t.m2,t.m3,t.m4,t.m5,t.m6)

## "Population (log)",
## "MID side B",      
## "Common law",      
## "GDP cap (log)",   
## "Polyarchy",       
## "MID side A"

t.combined$Variables.factor <- factor(t.combined$Variable,
                                      levels=c("Common law",
                                               "Population (log)",
                                               "GDP cap (log)",
                                               "Polyarchy",
                                               "Milex per cap",
                                               "MID"
                                               ),
                                      ordered=TRUE)

t.combined$type.factor <- factor(t.combined$type,
                                 levels=c("No force,\nfrom 1945",
                                          "Force,\nfrom 1945",
                                          "Overall,\nfrom 1945",
                                          "No force,\nfrom 1991",
                                          "Force,\nfrom 1991",
                                          "Overall,\nfrom 1991"),
                                 ordered=TRUE)

## http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")


## knitr FIGURE_MID_TOTAL
figure.a1 <- ggplot(t.combined, aes(color=type.factor)) +
    geom_hline(yintercept=0, color=gray(1/2), lty=2) +
#    geom_linerange(aes(x=Variables, ymin=X2.5., ymax=X97.5.),
#                       lwd=1, position=position_dodge(width = 1/2)) +
    geom_pointrange(aes(x=Variables.factor, y=b, ymin=X2.5., ymax=X97.5.),
                    lwd=2, position=position_dodge(width = 3/5),
                    shape=21, fill="WHITE") +
#    scale_color_brewer(palette="Set1") +
    scale_color_manual(values=cbbPalette) +
    scale_y_continuous(limits=c(-5,2),breaks=seq(-5,2,1)) + #name="SATT",
    coord_flip() +
    theme_bw() +
#    ggtitle("") +
    theme(axis.text.x = element_text(size = 25),#angle=45),
          axis.text.y = element_text(size = 25),
          axis.title.x = element_blank(),
#        axis.title.x = element_text("b", size=16),
          axis.title.y = element_blank(),
#          axis.title.y = element_text(size=16),
#          legend.title=element_blank(),
          legend.text=element_text(size=18),
          legend.key.size=unit(.5, "in"),
          legend.key.width=unit(1, "in"),
          legend.title=element_text(size = 18, face="bold"),
          legend.position="bottom",
          plot.title = element_text(lineheight=.8, face="bold",size=40)) +
    labs(color = "MID track record\n") +
    guides(color=guide_legend(nrow=2,byrow=TRUE,reverse=TRUE))

ggsave("gc-jg-project/kampala/paper/figures/app-figure-a1.png", plot = figure.a1, width = 14, height = 8, device = "png")


t.n.obs <- m1$n
t.n.fail <- m1$nevent
t.n.fail.tot <- table(Data$kampala)[2]
t.miss.vals <- dim(Data)[1]-t.n.obs

p.mid <- t.combined[t.combined$Variables=="MID","p"]
p.mid.below05 <- length(p.mid[p.mid<=.05])
p.mid.above05 <- length(p.mid[p.mid>.05])


p.val.overall.1991 <- round(m.m1[1,3],3)
p.val.overall.1945 <- round(m.m2[1,3],3)

p.val.force.1991 <- round(m.m3[1,3],3)
p.val.force.1945 <- round(m.m4[1,3],3)

p.val.no.force.1991 <- round(m.m5[1,3],3)
p.val.no.force.1945 <- round(m.m6[1,3],3)

p.wealth <- t.combined[t.combined$Variables=="GDP cap (log)","p"]
p.wealth.below05 <- length(p.wealth[p.wealth<=.05])
p.wealth.above05 <- length(p.wealth[p.wealth>.05])


ave.hazard.mid <- round(mean(c(exp(coef(m1)[1]),exp(coef(m2)[1]),exp(coef(m3)[1]),
                                     exp(coef(m4)[1]),exp(coef(m5)[1]),exp(coef(m6)[1]))),1)
ave.perc.hazard.mid <- 100*(1-ave.hazard.mid)


ave.hazard.milex <- round(mean(c(exp(coef(m1)[3]),exp(coef(m2)[3]),exp(coef(m3)[3]),
                                     exp(coef(m4)[3]),exp(coef(m5)[3]),exp(coef(m6)[3]))),1)
ave.perc.hazard.milex <- 100*(1-ave.hazard.milex)


ave.hazard.polyarchy <- round(mean(c(exp(coef(m1)[4]),exp(coef(m2)[4]),exp(coef(m3)[4]),
                                     exp(coef(m4)[4]),exp(coef(m5)[4]),exp(coef(m6)[4]))),1)
ave.perc.hazard.polyarchy <- 100*(ave.hazard.polyarchy-1)


ave.hazard.wealth <- round(mean(c(exp(coef(m1)[5]),exp(coef(m2)[5]),exp(coef(m3)[5]),
                                     exp(coef(m4)[5]),exp(coef(m5)[5]),exp(coef(m6)[5]))),1)
ave.perc.hazard.wealth <- 100*(ave.hazard.wealth-1)


ave.hazard.comlaw <- round(mean(c(exp(coef(m1)[7]),exp(coef(m2)[7]),exp(coef(m3)[7]),
                                     exp(coef(m4)[7]),exp(coef(m5)[7]),exp(coef(m6)[7]))),1)
ave.perc.hazard.comlaw <- 100*(1-ave.hazard.comlaw)


array_3d <- simplify2array(BF.mid.table)
mean_matrix <- apply(array_3d, 1:2, mean)
rownames(mean_matrix) <- paste0("Row", 1:6)
colnames(mean_matrix) <- c("H1", "H2", "H3")

BF.mid.table <- data.frame(mean_matrix)

BF.mid.table$varnames  <- c("Overall, from 1991",
                            "Overall, from 1945",
                            "Force, from 1991",
                            "Force, from 1945",
                            "No force, from 1991",
                            "No force, from 1945")

BF.mid.table <- BF.mid.table[c(1,3,5,2,4,6),c(4,1:3)]


colnames(BF.mid.table) <- c("MID track record",c("Pr($\\beta < 0$ | D)","Pr($\\beta = 0$ | D)","Pr($\\beta > 0$ | D)"))


BF.mid.table.print <- xtable(BF.mid.table,caption="Posterior Probabilities for the Coefficients in Figure~\\ref{app:overall-mid}",
                   label="app:posterior-overall-mid")
align(BF.mid.table.print) <- "rlccc"

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- c(-1) # toprule
addtorow$pos[[2]] <- c(0) # midrule
addtorow$pos[[3]] <- c(nrow(BF.mid.table)) # endrule

addtorow$command <- c('\\toprule\n'
                     ,'\\midrule\n'
                     ,'\\bottomrule\n')



## @knitr BF_MID_TABLE_PRINT
print(BF.mid.table.print, sanitize.colnames.function=identity,
      sanitize.text.function = function(x) {x},
      floating = TRUE, #floating.environment = "sidewaystable",
#      scalebox=.87,
      booktabs = TRUE,
      table.placement = "ht!",
      caption.placement = "top",
      add.to.row = addtorow,
      include.rownames=FALSE,
      hline.after=NULL
#      ,size="\\scriptsize"
#     ,file="regression/paper/tables/tab-reg-models.tex"
      )
